#include <lversion.h>
#if L_VERSION_MODULE_EQUATION == 111225

#include "_lequation_111225.h"

#include <boost/assert.hpp>

using namespace llib;
using namespace llib::math;
using namespace std;

// cubic equation solution id
// 0: the discriminant with complex number
// 1: the discriminant with trigonometric function

/* Log:
 * [2011-11-24]: switch to using Double floating variables for internal precision.
 */

/* NOTE:When solving the cubic equation with three real roots,  */
/*      solution 1 is about THREE times faster than solution 0! */
/*      A simple experiment is done at 20110623                 */
#define CUBIC_EQUATION_SOLUTION	1

#if CUBIC_EQUATION_SOLUTION == 0
#include <complex>
#endif

vector<Float> EquationSolver::SolveCubicEquation_base(const Float A, const Float B, const Float C, const Float D)
{
    Double b = B / A;
    Double c = C / A;
    Double d = D / A;
    // local variable for notation & compute convenience
    // q = c/(3*a) - b^2/(9*a^2)
    // r = (b*c)/(6*a^2) - b^3/(27*a^3) - d/(2*a)
    Double q = (3 * c - b * b) / 9;
    Double r = -27 * d + b * (9 * c - 2 * b * b);
    r /= 54;
    // discriminant = q^3 + r^2
    Double disc = q * q * q + r * r;
    // the left part
    Double term = b / 3;

    if(abs(disc) < EPS_F) {
        // when the discriminant is equal to zero, there's at least two equal roots
        if(abs(r) < EPS_F) {
            // if r is equal to zero, there's three equal roots
            return vector<Float>(1, Float(-term));
        } else {
            // a single root with two equal roots
            vector<Float> res(2);
            if(r > 0) {
                Double r1_3 = pow(r, 1.0/3.0);
                res[0] = Float(-r1_3 - term);
                res[1] = Float(r1_3 * 2 - term);
            } else {
                Double r1_3 = -pow(-r,1.0/3.0);
                res[0] = Float(r1_3 * 2 - term);
                res[1] = Float(-r1_3 - term);
            }
            return res;
        }
    } else if(disc > 0) {
        // when the discriminant is greater than zero, just one real root
        Double discSqrt = sqrt(disc);
        Double s = r + discSqrt;
        s = (s < 0 ? -pow(-s, 1.0/3.0) : pow(s, 1.0/3.0));
        Double t = r - discSqrt;
        t = (t < 0 ? -pow(-t, 1.0/3.0) : pow(t, 1.0/3.0));
        return vector<Float>(1, Float(s + t - term));
    } else {
        // when the discriminant is less than zero, all roots are real and unequal
        Float x0, x1, x2;
#if CUBIC_EQUATION_SOLUTION == 0
        // here we introduce complex number
        complex<Double> discComplex(disc);
        complex<Double> discSqrtComplex  = sqrt(discComplex);
        complex<Double> leftComplex      = pow(r + discSqrtComplex, 1.0/3.0);
        complex<Double> rightComplex     = pow(r - discSqrtComplex, 1.0/3.0);
        complex<Double> w0(-0.5f,  SQRT3/2);
        complex<Double> w1(-0.5f, -SQRT3/2);
        // x0
        complex<Double> tmpComplex = leftComplex + rightComplex;
        BOOST_ASSERT( abs(tmpComplex.imag()) < EPS_F );
        x0 = Float(tmpComplex.real() - term);
        // x1
        tmpComplex = w0 * leftComplex + w1 * rightComplex;
        BOOST_ASSERT( abs(tmpComplex.imag()) < EPS_F );
        x1 = Float(tmpComplex.real() - term);
        // x2
        tmpComplex = w1 * leftComplex + w0 * rightComplex;
        BOOST_ASSERT( abs(tmpComplex.imag()) < EPS_F );
        x2 = Float(tmpComplex.real() - term);
#else if CUBIC_EQUATION_SOLUTION == 1
        // use trigonometric functions
        /* 20110624:Experiment shown that the functions:    */
        /*  (WHY?)  acos, sqrt, cos is about 3% faster than */
        /*          the float version: acosf, sqrtf, cosf.  */
        q = -q;
        Double q_3 = q * q * q;
        Double angle = r / sqrt(q_3);
        BOOST_ASSERT( angle >= -1 && angle <= 1 );
        Double dum = acos(angle);
        Double r13 = 2.0f * sqrt(q);
        x0 = Float(r13 * cos(dum / 3.0f)              - term);
        x1 = Float(r13 * cos((dum + 2.0 * PI) / 3.0f) - term);
        x2 = Float(r13 * cos((dum + 4.0 * PI) / 3.0f) - term);
#endif // CUBIC_EQUATION_SOLUTION
        // result sort
        vector<Float> res(3);
        if(x0 < x1) {
            if(x0 < x2) {
                if(x1 < x2) {
                    res[0] = x0; res[1] = x1; res[2] = x2;
                } else {
                    res[0] = x0; res[1] = x2; res[2] = x1;
                }
            } else {
                res[0] = x2; res[1] = x0; res[2] = x1;
            }
        } else {
            if(x0 < x2) {
                res[0] = x1; res[1] = x0; res[2] = x2;
            } else {
                if(x1 < x2) {
                    res[0] = x1; res[1] = x2; res[2] = x0;
                } else {
                    res[0] = x2; res[1] = x1; res[2] = x0;
                }
            }
        }
        return res;
    }
}

vector<Float> EquationSolver::SolveQuadraticEquation_base(const Float A, const Float B, const Float C)
{
    // the discriminant
    Double delta = B * B - 4 * A * C;
    if(abs(delta) < EPS_F) {
        return vector<Float>(1, -B / A / 2);
    } else if(delta > 0) {
        Double deltaSqrt = sqrt(delta);
        vector<Float> res;
        res.push_back(Float((-B - deltaSqrt) / A / 2));
        res.push_back(Float((-B + deltaSqrt) / A / 2));
        if(A < 0) swap(res[0],res[1]); // [2012-1-9] : we want to sort the result, but left the bug till now
        return res;
    } else {
        return vector<Float>();
    }
}

vector<Float> EquationSolver::SolveLinearEquation_base(const Float A, const Float B)
{
    return vector<Float>(1, -B / A);
}

vector<Float> EquationSolver::SolveCubicEquation( const Float A, const Float B, const Float C, const Float D )
{
    if(abs(A) < EPS_F) {
        // degenerate to quadratic: B*x^2 + C*x + D = 0
        return SolveQuadraticEquation(B, C, D);
    } else if(abs(D) < EPS_F) {
        //  factorization: x * (A*x^2 + B*x + C) = 0
        vector<Float> res = SolveQuadraticEquation(A, B, C);
        Int32 I = 0;
        for(; I < SZ(res); ++I) {
            if(abs(res[I]) < EPS_F) break;
        }
        if(I == SZ(res)) {
            // res has no zero value inside
            res.push_back(0.f);
        }
        return res;
    } else {
        return SolveCubicEquation_base(A, B, C, D);
    }
}

vector<Float> EquationSolver::SolveQuadraticEquation( const Float A, const Float B, const Float C )
{
    if(abs(A) < EPS_F) {
        // degenerate to linear: B*x + C = 0
        return SolveLinearEquation(B, C);
    } else if(abs(C) < EPS_F) {
        // factorization: x * (A*x + B) = 0
        vector<Float> res = SolveLinearEquation(A, B);
        if(SZ(res) > 0 && abs(res[0]) > EPS_F) {
            res.push_back(0.f);
        }
        return res;
    } else {
        return SolveQuadraticEquation_base(A, B, C);
    }
}

vector<Float> EquationSolver::SolveLinearEquation( const Float A, const Float B )
{
    if(abs(A) < EPS_F) {
        // return an empty result
        // we treat both no root and infinite root as invalid
        return vector<Float>();
    } else {
        return SolveLinearEquation_base(A, B);
    }
}

#endif
